Superdiffusive Heat Transport in a class of 
Deterministic One-Dimensional Many-Particle 

Lorentz gases 



Pierre Collet 1 , Jean-Pierre Eckmann 23 and 
Carlos Mejia-Monasterio 4 

Centre de Physique Theorique, CNRS UMR 7644, Ecole Poly- 
technique, F-91128 Palaiseau Cedex (France) 
2 Departement de Physique Theorique, Universite de Geneve 
3 Section de Mathematiques, Universite de Geneve 
4 Istituto dei Sistemi Complessi, Consiglio Nazionale delle 
Ricerche, Sesto Fiorentino Italy 

Abstract 

We study heat transport in a one-dimensional chain of a finite number N of 
identical cells, coupled at its boundaries to stochastic particle reservoirs. At the 
center of each cell, tracer particles collide with fixed scatterers, exchanging mo- 
mentum. In a recent paper, [1], a spatially continuous version of this model 
was derived in a scaling regime where the scattering probability of the tracers 
is 7 ~ 1/N, corresponding to the Grad limit. A Boltzmann type equation de- 
scribing the transport of heat was obtained. In this paper, we show numerically 
that the Boltzmann description obtained in [1] is indeed a bona fide limit of the 
particle model. Furthermore, we also study the heat transport of the model when 
the scattering probability is one, corresponding to deterministic dynamics. At 
a coarse grained level the model behaves as a persistent random walker with a 
broad waiting time distribution and strong correlations associated to the deter- 
ministic scattering. We show, that, in spite of the absence of global conserved 
quantities, the model leads to a superdiffusive heat transport. 
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1 Introduction 

The rigorous description of (classical) non-equilibrium steady states (NESS) re- 
mains an elusive problem, despite some remarkable progress in the last few years, 
as described, e.g., in [2, 3]. One of the main reasons seems to be the impossibility 
to guess the stationary state, which is one of the mechanisms which work in equi- 
librium statistical mechanics. There are therefore many studies which try to un- 
derstand better what the essentials of NESS are, how different models fit together, 
and what are the best descriptions of NESS. Among the few works where the 
NESS have been obtained, we mention the harmonic chain coupled to Langevin 
heat baths [4] and anharmonic chains coupled to infinite (non compact) reservoirs 
[5,6]. 

In this paper, we study heat and particle conduction in a 1 -dimensional model 
introduced in [1], and which is a variant of a model studied in [7, 8]. We ana- 
lyze several of its properties and in particular show, by numerical study, that the 
Boltzmann description of the model obtained in [1] is indeed a bona fide limit 
of the particle models studied here, but only so if the coupling, per site, is of or- 
der 0(1/ N) when the number of sites is (large) N. The model does not seem to 
obey the Fourier law. The reader familiar with transport problems will also notice 
that (due to the strictly 1 -dimensional character of the model), the particle num- 
ber in the NESS must be infinite. However, most of the particles will have very 
small kinetic energy, so that the system has very nice energy profiles, which, so to 
speak, are generated by those particles with energy away from 0. This is in fact 
reminiscent of non-normalizable measures in dynamical systems [9]. 

The model in question consists of a chain of identical cells, each of which con- 
tains a fixed point-like scatterer that exchanges momentum with tracer particles. 
Inside the system the particles move deterministically between cells, interacting 
with the scatterers but not among themselves. However, on their passage the par- 
ticles modify the local state of the substrate, which in turns alters the evolution of 
the other particles. At the collisions energy is conserved. At its boundaries the 
chain is in contact with two stochastic particle reservoirs, characterized by a fixed 
temperature. 

In the next section we describe in detail the model and in Sect. 3 the general 
properties of its non-equilibrium steady state. In Sect. 4 we study the continuous 
limit of our model and compare it with the theory that appeared in [1]. Finally, 
in Sects. 5 and 6 we discuss the energy and particle transport of the deterministic 
finite chain. 
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2 The Model 



In this section, we describe the 1 -dimensional particle model that we consider, 
and briefly review the results of [1] for the continuous model. 

The model consists of N cells in a row, each cell of length A. In the center of 
each cell there is a point-like scatterer which does not move but which has a "mo- 
mentum" P 6 1 and a mass M. Particles move in these cells. They have mass 
m 7^ M and momentum p. The particles do not interact among themselves but 
they do interact with the scatterers as follows: Whenever a particle with momen- 
tum p reaches a scatterer whose momentum is P, the following happens: With 
probability 1 — 7/iV, the particle crosses to the other side of the scatterer, and 
continues with momentum p, while the scatterer retains its momentum P. With 
probability 7/iV actual scattering takes place and the new momenta p and P are 
given by 



When 7 = N, the particles interact with the scatterer every time they encounter 
one, and the model is fully deterministic, except for the nature of the baths. These 
rules are similar in spirit (but far more rich), to the flipping Lorentz lattice gases 
studied some years ago [10]. When 7 < N, then the model has some randomness, 
since we need to decide whether scattering takes place, or the particle flies through 
the scatterer. Note, however, that this randomness does not change the energies of 
the actors in the system. 

The collision rules are just those of elastic scattering, but with the scatterers 
not moving. The deterministic model (with 7 = N) is a one dimensional general- 
ization of the models previously studied in [7, 8, 11, 12], where the scatterers are 
fixed freely rotating disks. In these models, the collisions provide a local energy 
mixing among different degrees of freedom that leads to a local state which is well 
approximated by a local equilibrium state where the particles behave as a perfect 
gas. Close to equilibrium, Green-Kubo relations for the heat and particle fluxes 
are valid and the corresponding Onsager reciprocity relations are satisfied [7, 11]. 
Moreover, in the zero-coupling limit, where the invariant measure of the NESS 
is expected to be multivariate Gaussian, compact analytical expressions for the 




where the scattering matrix S is 




(2.1) 
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currents and the density profiles have been obtained [8]. To the lowest order, the 
corrections due to a finite coupling were considered in [12]. Note that if particles 
and scatterers have the same mass, the model has an infinite number of conserved 
quantities, as their momenta are just exchanged in the collisions. Trajectories for 
one single particle have been considered in [13]. 

To force the system out of equilibrium, we couple the leftmost and rightmost 
cells to infinite ideal particle reservoirs. From each reservoir particles are injected 
into the system at a given rate v and with momenta distributed according to 

F(p)dp = e(±p)(27rmk B T)- 1/2 e~ p2 / 2mk * T dp , (2.2) 

where T is the temperature of the reservoir, k B the Boltzmann constant and the 
Heaviside function 0(p) restricts the sign of the momentum according to the side 
from where the particles are injected ("— " for those entering from the right side 
and "+" for those entering from the left). Equation (2.2) implies that the particles 
in the vicinity of the opening between the system and the reservoir have a momen- 
tum density f(p) with a non-normalizable singularity at p = 0. For the reservoir 
coupled to the leftmost (i = 1) cell, f(p) is 

kip) dp = — , n dp , (2.3) 

(2irmk B T L ) l/2 \p\ 

and for the reservoir coupled to the rightmost (i = N) cell 

0(—v) e ~p 2 /2mk B T R 

/r(p) dp = - ;\ 1/2 n dp . (2.4) 

(2nmk B T R ) L/2 \p\ 

The reader should note that the singularity of f(p) is non-integrable only for 
one dimensional models, more precisely, for models with one dimensional dy- 
namics. This is related to the fact that, while in any dimension, slow particles 
need much longer times to move across the system than fast particles, only in 
1 -dimensional dynamics does this imply that particles steadily accumulate near 
p — 0. In higher dimensions, l/\p\ is integrable near p = 0. Thus, the particle 
density diverges in the stationary state as t — > oo. Due to this seemingly un- 
physical property, in the past, it has been argued that injecting particles with the 
momentum distribution (2.2) cannot be correct [14]. However, as has been shown 
in [1], it is the reservoir distribution f(p) and not F(p) that admits stationary solu- 
tions which, at the same time, preserve the distribution of the scatterers' momenta. 
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Furthermore, the higher momenta of p are obviously well-behaved in the limit of 

t — > oo. 

More precisely, denoting by g(P^) the distribution of the momentum of the i-th. 
scatterer, it was shown in [1] that at equilibrium (T L = T R = T), both F(p) and g 
are Gaussian and are given by 

F(p)dp = {2itmk B T)- l l 2 e-P l 2mk * T dp , (2.5) 

g (P) dP = (2iTMk B T)- 1 / 2 e - p2 /2Mfc B T dp _ (2.6) 

As mentioned before, because of the factor l/\p\ in /, the number of particles 
in the system is infinite at stationarity. Starting with any initial distribution, and 
letting the system evolve in time t, the number of particles with low momenta 
grows without bounds. On the other hand if one can prove that the number of 
particles with momentum p > p has a limit as t — > oo then the stationary state 
is well defined, in spite of the divergence of the total number of particles [9]. To 
convince ourselves that this is indeed the case, we have measured the evolution 
in time of the number of particles in the system and their momentum distribution, 
for a chain of N = 201 cells at equilibrium: particles were injected from both 
reservoirs at the same rate v = 100 with the same temperature T = 100. The 
results are shown in Fig. 1, where the number of particles n t {p) with momentum 
p times the momentum \p\, is plotted for several successive times. We observe that 
for p larger than some p (t), pn t (p) is stationary. This means that the total number 
of particles nit) = n t (p) dp diverges due to the slow accumulation of cold 
particles. Indeed, n(t) diverges logarithmically with time, as can be seen in the 
inset of Fig. 1. In the same vein, p (t) is seen to decrease to zero logarithmically. 

In the limit of N — > oo, setting x = i/(N ■ X) for the position of the i-th 
cell, g(Pi) — > g(P, x), is the probability density that the scatterer at position x has 
momentum P and so J dP g(P, x) = 1, by definition. In [1] it was argued that in 
the continuous limit, this system can be modeled by a Boltzmann equation, whose 
stationary solution is described by the equations 

pd x F(p,x) = 'y\p\ J dP (F(p, x)g(P , x) - F(p, x)g(P , x)j , (2.7a) 
= J dp (F(p, x)g(P , x) - F(p, x)g(P , x)) , (2.7b) 

where the quantity F(p, x) is equal to F(p, x) = \p\f(p, x), with f(p, x) the prob- 
ability density that a particle at position x has momentum p, and 7 e [0, JV] corre- 
sponds to the same quantity of the particle model, which determines the scattering 
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Figure 1: Number of particles n t (p) with momentum « p times the momentum p, 
at t = 200, 600, 1000, 1400 and 1800, for a chain of N = 201 cells at equilibrium 
with 2/ L = u R = 100 and T L = T R = 100, a = 1/2 and 7 = N. dp is the width 
of the bins used to compute the empirical distribution n t (p). Thus, the s-axis 
corresponds to the bin number. In the inset: logarithmic divergence of the total 
number of particles inside the system. 



probability. Therefore, for p > 0, F(p, x) is related to the rate of particles with 
momentum p moving from x to the right. Similarly, F(p, x) is, for p < 0, related 
to the rate of particles with momentum p moving from x to the left. Note that, in 
contrast to fip), the function F(p) is free of singularities. Therefore, in spite of 
the infinite number of particles in the stationary state of the scattering model, the 
flux F(p, x) is finite and integrable, and it is for this quantity that the Boltzmann 
equation is formulated. 

The Boltzmann equation (2.7) was derived assuming that F(p) and g(P) are 
statistically independent. Therefore, the similarity between the particle model and 
its Boltzmann version should be best in the case of large iV and when there are 
many particles (with momentum \p\ > p > 0) in each cell. Moreover, in [1], 
it was proven that for any particle injections / L (p) and f R (p) in a certain cone in 
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Banach space, (2.7) has solutions when 7 = (9(1). Our numerical studies are, 
however, for parameter values well outside this cone, and still give a very good 
comparison between the particle model and the Boltzmann model. Furthermore, 
when 7 = O(N), the particle model is still well defined, although different from 
the Boltzmann model. If 7 = N a particle will scatter whenever it meets a scat- 
terer and will never fly to the other side of a scatterer without collision. When 
7 < IV, the local state cannot correspond to a local equilibrium state. Indeed, 
when scattering is rare, the particles do not interact, thus leading to local states 
that are described by the sum of two different families of particles: those that 
were injected from the left, flying ballistically to the right, and those injected 
from the right that fly ballistically to the left [15]. However, if 7 = N, then all 
particles scatter when they encounter a scatterer, leading to a stronger local in- 
teraction. One would expect that this strong interaction would lead to a diffusive 
particle behavior. However, in Sect. 5, we will show that the transport remains 
superdiffusive. 

3 Non-equilibrium steady state 

In this section we consider the model described in the previous section, coupled 
to reservoirs injecting particles into the chain, at different rates and at different 
temperatures. In the rest of the paper we will study our model for 7 = 1 (corre- 
sponding to the Grad limit considered in [1]) and for 7 = iV (corresponding to the 
deterministic particle dynamics). 

Out of equilibrium, particle and energy currents appear, whose magnitudes are 
determined by the differences of injection rates and of temperatures. The particle 
injection rate is defined as 



where F(p) is the momentum distribution of the injected particles, given by (2.2). 
As we have discussed, the particle density in the bath is infinite. However, since 
the integral in (3.1) is finite, one realizes that v takes the role of an effective 
particle density. Moreover, since the injected particles correspond to a perfect gas 
at equilibrium, one can define the chemical potential of the reservoir as 




(3.1) 




(3.2) 



with /jl a pure constant. The injection rate v also trivially determines the rate at 
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Figure 2: Profiles of the particle density (lower panel) and the kinetic energy per 
particle density E K (upper panel), for a chain of N = 30 cells, a = 0.5 and for 
7 = 1 (squares), and 7 = N (circles). Particles were injected to the chain with 
rates z/ L = 220, u R = 180 and temperatures T L = 81.818 and T R = 122.222, 
indicated here by the dashed line. Here, T(x) is computed as the mean scatterers 
kinetic. 



which energy is injected into the system as 

e = uk B T . (3.3) 

We now proceed to study the non-equilibrium state of our model. In Fig. 2 
we compare the particle density g(x) and the temperature profiles numerically 
obtained for two different values of the scattering parameter, 7 = 1 and 7 = N. 
The temperature at the i-th cell is computed as the average p 2 with respect to 
the particle momentum distribution function Fj(p) measured at the i-th cell. This 
temperature coincides with the time averaged kinetic energy of the i-th scatterer, 
indicating a good local equilibration. All the observables were averaged for a 
time interval during which the total number of particles inside the channel does 
not change appreciably. 
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The temperature in the bulk of the system does not match the nominal temper- 
atures of the reservoirs (indicated by the dashed lines in the upper panel). More- 
over, we observe that the energy mismatch at the contact with the reservoirs de- 
pends on the scattering probability. For 7 = 1 the temperature profile is less steep 
than for 7 = N as expected, since in the former case, particle-particle interaction, 
mediated by the scatterers, is less effective than in the latter case. As for the par- 
ticle density, a smaller scattering probability leads to a less steep density profile. 
Moreover, the accumulation of particles at the contact with the reservoirs is more 
pronounced when the scattering is less frequent. This is because the scattering 
contributes to heat up the injected cold particles. In the limit of iV — > 00, the 
accumulation of particles at the boundaries produces singularities [14]. 

We also have measured the rate at which particles cross from one cell to an- 
other from left to right v R (x) and from right to left f L (x). 

Out of equilibrium, the particles inside the cell at position x, leave the cell 
to the right at a different rate than the rate at which they leave the cell to the 
left. This is a consequence of the substrate-mediated particle-particle interaction. 
Calling the rate at which particles cross from one cell to another from left to right 
v R (x) and from right to left u L (x), the particle current is defined as 

J n (x) = v R (x) - v L (x + dx) . (3.4) 

Therefore, these local rates are strongly dynamically constrained, so that the sta- 
tionary particle current is uniform. Indeed, we find that, in the stationary state, 
a uniform particle current, with extremely linear profiles for u R (x) and u L (x). In 
Sect. 5, we will see that this strong correlations determine an unexpected superdif- 
fusive particle and energy transport. 

4 Infinite volume limit: 7=1 

In this section, we compare the non-equilibrium probability distribution functions 
of the discrete model with those predicted by the Boltzmann equation (2.7). 

We have solved numerically (2.7) by a discretization in momentum space. 
Fixing the spacing of the discretized momentum to Ap, the number of points 
we considered is the minimum necessary to keep the information from the tails 
of the distributions as small as ~ 10~ 10 . We proceed as follows: at any x, the 
equation (2.7b) for g only depends on F(-, x). We discretize (2.7b) and solve it as 
an eigenvalue problem (with eigenvector g(-, x)). The only limitation is the size 
of the matrices one obtains in this way: Our runs were done with matrices of size 
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Figure 3: Momentum distribution F(p, x) of the particles at the left (solid curve) 
and right (dashed curve) ends of the system , from the solution of the Boltzmann 
equation (2.7), for a = 0.5, u L = 220, z/ R = 180, T L = 81.818 and T R = 122.222. 
In the inset the finite size deviation of F N (p, x) from the Boltzmann solution are 
shown for iV = 2 (black), N = 4 (red), N = 8 (green) and N = 16 (blue). 



~ 4200. The function g found in this way is then inserted into (2.7a). High-order 
integration in position space is then used to integrate (2.7a) from x = to x = 1. 
Therefore, fixing the injection at x = to (2.3), we use a shooting method to 
determine the extraction of particles at x = in such a way that at x = 1 the 
desired injection (2.4) will result (see also [1] for more details). 

In Fig. 3 we show the solution of (2.7) for a = 0.5, v L = 220, v R = 180, 
T L = 81.818 and T R = 122.222, and 7 = 1. The first peculiarity of the non- 
equilibrium distributions is the jump at p — 0. This is due to the very weak 
particle-particle interaction obtained for 7 = 1. The size of the jump is partly 
determined by 7 and, as it is clear from Eq. (3.1), partly by the difference between 
u L and z/ R . Note that, as a consequence of the temperature gradient, the positive 
and negative parts of the distribution are only approximately Gaussian. They are 
Gaussian only if T L = T R . 
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Figure 4: Momentum distribution F N (p,x) of the particles at x = l/N (solid 
curve) and x — 1 (dashed curve), for a chain of N = 30 cells, with o = 0.5, for 
a) 7 = 1 and b) 7 = N. The bath's parameters are as in Fig. 3. 



To study the limit N — > 00, we have numerically followed the evolution of 
finite size iV chains and measured the particle momentum distribution F N (p, x) 
for the same parameters as above. In the inset of Fig. 3, the deviation of F N (j>, x) 
from the solution of the Boltzmann equation F(p, x) is shown for chains from 
N = 2 to iV = 16. 

Deviations are seen over the whole domain, although they are biggest at the 
center of the distribution. They are not symmetric in p, which is an indication that, 
for a given size N, deviations may depend on the injection rates and bath's tem- 
peratures in general. Furthermore, we observe that the solution of (2.7) appears 
to be the asymptotic distribution, lim^^^ F N (p, x) — > F(p, x). In any case, the 
deviations from F(p, x) are less than 0.1%, tending to zero very fast. For instance, 
for a chain of N = 16 the deviations are less than 0.01%. 

Finally, we have also measured the distribution F N (p, x) for the deterministic 
finite chain (7 = N). In Fig. 4, we show F N (j>, x) at x = l/N (solid curve) and 
x = 1 (dashed curve), for a chain of iV = 30 cells and 7 = 1 (panel a), 7 = N 
(panel b). The other parameters are reported in the caption of Fig. 3. The distribu- 
tions F 30 (p, x) in Fig. 4-a, are on top of the solution (2.7) and, as mentioned above, 
the deviations from the asymptotic distribution decay very fast. For 7 = N, the 
jump at p = is much smaller, with the only remaining the contribution coming 
from the difference |z/ L — v R \ of the injection rates. As expected, the distributions 
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in Fig. 4-6 are not Gaussian for this non-equilibrium case. 
5 Energy transport 

In this section we turn our attention to the heat transport of the deterministic 
model, i.e., setting 7 = N. The particles of mass m (which we call tracers) 
are the only energy carriers of the system. We start by analyzing their dynamics. 

5.1 Microscopic evolution 

We consider a finite chain of N cells with periodic boundary conditions and n 
particles per cell. A stationary state of this closed system is the equilibrium state 
characterized by N, n and the total energy E given by 



We start the evolution with the scatterers at momenta. As for the open chain, 
the state of the system approaches equilibrium logarithmically slowly (in time). 
During the transient, the substrate continuously extracts energy from the gas of 
particles; the total energy of the scatterers grows logarithmically in time, until it 
saturates at sufficiently long times. All measurements are taken after the system 
has relaxed to the approximate equilibrium state. 

Once equilibrium is reached, we focus on the evolution of a tagged particle in 
the bulk of the system, when N and n are sufficiently large. This is convenient, 
since particles interact among themselves only through their collisions with the 
substrate, and thus the local dynamics depends on the particle density. Here, we 
are interested in the high density regime, which, following the discussion about 
Fig. 1, is a good approximation of the stationary n = 00 state. 

When the particle encounters a scatterer, its velocity after collision is deter- 
mined by (2.1). In fact, this scattering matrix leads to a persistent motion of the 
particle, namely the probability that the particle continues in the same direction in 
which it reached the scatterer, is larger than 1/2. This probability, fi, can be easily 
computed as follows 1 : without loss of generality assume that before the collision, 
the particle's velocity is v > 0. In equilibrium, the scatterer velocity distribution 
function is obtained by taking V = P/M in (2.6). Taking into account that af- 
ter the collision the particle's velocity is v' = —ov + (1 + o)V (see (2.1)), the 

'To analyze the dependence of fi on the masses, it is convenient to work with the velocities 
instead of the momenta. 




(5.1) 
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probability that after the collision the particle has a velocity v' > can be written 
as 

(mM) 1/2 r°° rnv 2 f 00 MV 2 

H(a) = P(v' > 0\v > 0) = - — — / dv e / dVe ™*r , (5.2) 

2vr£; B T J J * v 

that can be integrated to yield 

/x 1 / m \ 1/2 f°° , /(M-m>\ -^i 

where erf(-) is the error function. The limit values of /i can be easily read from 
(5.3): for a = — 1, namely M — 0, the error function is erf(— oo) = — 1 and 
li — 1. In the opposite case, when a = 1 (M = oo), /i = 0. Finally, for cr = 0, 
namely M = m, erf(0) = and fj, — 1/2. With the exception of er = 0, the 
dynamics of the particles is persistent. In Fig. 5, the probability fi(a), computed 
from the statistics of the collisions of the tagged particle is shown. 

At a coarse grained description, the dynamics of the particles can be seen as a 
persistent random walk with waiting time r corresponding to the collision times, 
that are determined by the particle's velocity. We have measured the distribution 
of the waiting time *(r) of the tagged particle for different values of a. In Fig. 6 
we show \&(t) for a = and a = 0.5. As a consequence of the single particle's 
velocity distribution, ^(r) turns out to be a broad distribution \l/(r) ~ t~ (1+s) , 
with s ~ 1 for a = 0.5 and s ~ 2 for a = 0. Therefore, our persistent walker 
seemingly performs a Levy walk. The power —2 for a ^ can be derived (ap- 
proximately) from a multiple integral as in (5.3). 

In the continuous limit, a persistent random walker yields to a particle's den- 
sity whose evolution is described by the telegraph equation [16]. Noting that 
asymptotically the telegraph equation yields to a diffusive evolution and that the 
Levy walk for 1 < s < 2 does only induce anomalous corrections to the normal 
long-time behavior [17], one would expect that the microscopic particle's dynam- 
ics yields diffusive transport. However, this is not the case. In Fig. 7 we show the 
evolution of the dispersion of the position of a tagged particle (x 2 (t)} for a = 0.5, 
averaged over an ensemble of initial conditions. Asymptotically, (x 2 (t)) ~ t a , 
with a < 2. In fact, we have found that the asymptotic scaling of (x 2 (t)) depends 
on the mass ratio parameter a (see inset of Fig. 7). For a = the particle's motion 



2 There is no factor l/\v\ here, because we must consider the probability of a particle with 
velocity in [v, v + dv] hitting a scatterer within a given time. 
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Figure 5: Persistent probability /i as a function of the mass ratio parameter a, 
averaged over the evolution of a tagged particle in a chain of N = 11, n = 20 and 
T = 500. 



is ballistic, while for <j / 0, the motion is superdiffusive. The observed anoma- 
lous behavior proves that in one dimension, the effect of the dynamical memory 
of the deterministic model is much stronger than in higher dimensions 3 . Since 
the particles are the energy carriers one expects that the energy transport will be 
anomalous as well. We study this in the next section. 

6 Heat conductivity 

We turn our attention to the energy transport of our model. Considering the open 
system coupled at its boundaries to two particle reservoirs, we have computed 
the dependence of the heat conductivity k on the size of the system N, for fixed 
nominal values of the injections and temperatures of the particle reservoirs. We 

3 As a note aside, if the direction of the particle after the collision if chosen randomly so that 
the effects of the dynamical memory can be neglected, then the diffusive transport is recovered. 
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Figure 6: Distribution function of the collision times ^(r), for a = and a = 0.5. 
The dashed lines correspond to fit to a power law. We obtain, in within numerical 
accuracy, ^(r) ~ t -3 for a = and ^(r) ~ r~ 2 for a = 0.5. 



define the heat conductivity as 



where J v is the measured energy current and T x (resp. T N ) is the temperature 
measured in the leftmost (resp. rightmost) cell that, as we have seen, in general 
does not coincide with the temperatures of the reservoirs. (The length of the 
system is 1 since we space the scatterers by A = 1/N.) The results of simulations 
are shown in Fig. 8 for 7 = 1 (squares) and 7 = N (circles). When the scattering 
is rare (7 = 1), we obtain k ~ N. This can be understood from the fact that 
particles move ballistically, interacting with the lattice very rarely. 

Surprisingly, the energy transport in the deterministic model (7 = N) is 
anomalous, with a heat conductivity that diverges as k ~ iV 1 / 3 . It is interest- 
ing to note that usually anomalous heat conduction is related to the existence of 
additional global conserved quantities [3]. However, for cr/0 the bulk dynamics 
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Figure 7: Time dependence of the variance of the position of the tagged particle 
(x 2 (t)), for N = 11, n = 20 and temperature T = 500 (solid curve), when a = 1. 
The dashed line corresponds to a scaling ~ t 2 . In the inset, asymptotic power 
scaling (x 2 (t)) ~ t a , as a function of the mass ratio parameter a. 



of our model only conserves energy. As far as we know, this is the first exam- 
ple of a mechanical model that, not having additional integrals of motion, shows 
anomalous heat conduction. 

6.1 Return to equilibrium 

In order to shed more light on the anomalous heat transport we studied the sys- 
tem's equilibrium response to a finite energy perturbation. Suppose that at a cer- 
tain initial time, t = 0, the equilibrium state of the system is perturbed by an 
additional amount of energy AE that is distributed among all the degrees of free- 
dom in a finite region of volume V, around the position x. By measuring the 
evolution of the energy field, one can estimate how heat propagates through the 
system. 

Considering the closed system as in Sect. 5.1 we proceed as follows: at time 
t = we perturb the state 5° of the system to S, as follows: the energies of the 
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Figure 8: Heat conductivity k, as a function of the number of cells N, for two 
values of 7: 1 (squares) and N (circles). The rest of the system's parameters are 
as in Fig. 3. For 7 = N, k diverges as ~ iV 1 / 3 , while for 7 = 1, k diverges as 
~ N. These power laws are indicated by the dashed lines. 



particles and scatterers contained in the J\f central cells are changed so that the 
total energy inside these cells is E pen . After this, we let the central subsystem 
relax. To obtain the evolution of the energy perturbation AE(x, t), we have fol- 
lowed two trajectories of the system: the unperturbed one, with initial state S° 
and the perturbed one, with initial state S. Then, the energy difference at time t 
and position x is 

AE(x = iX, t) = (E^t) - Ef(t)) , (6.1) 

where E^t) is the energy contained in the ith cell at time t of the perturbed trajec- 
tory and respectively for Ef(t), and (•) denotes the average over an ensemble of 
different initial realizations. 

When dynamical correlations are not too strong, one expects that after a suf- 
ficiently long time the perturbation AE(x, t) scales (with x measured from the 
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Figure 9: Evolution of the energy difference AE(x,t) for a chain of N = 101, 
n = 50 and a = 0.5, and for 7 = iV (left panel) and 7 = 1 (right panel). The 
initial energy per degree of freedom was e° = 5000 and s = 50000. 



initially perturbed cells) as 

AE(x,t) = ^AE (^,t) , (6.2) 

where the power £ is related to the scaling of the heat conductivity with the size 
of the system L as [18]: 

K = JV2-VC . (6.3) 

In particular, £ = 1/2 corresponds to normal diffusion, while £ = 1 corresponds 
to ballistic motion. 

In Fig. 9 we show the evolution of the energy difference AE(x, t) for the 
deterministic chain 7 = N (left panel) and compare it with the stochastic chain 
with 7=1 (right panel). We observe that for 7 = N, the initial excess of energy at 
the central cell decays very rapidly. Practically none of the initial energy remains 
in the central cell. The perturbation moves symmetrically to the ends of the chain, 
carried by two seemingly independent families of particles, those with positive 
velocity and those with negative velocity. Actually, it is the fast decay of the 
energy at the center that marks the existence of very strong dynamical correlations. 
A similar observation has been made recently in a random walk with memory 
in the waiting times of successive steps [19]. For 7 = 1, we also observe an 
initial fast decrease of the energy perturbation. The evolution of the peaks with 
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t t 

Figure 10: Damping of the maximum value of AE(x, t) as a function of time, for 
the same simulation as in Fig. 9, for 7 = N (left panel) and 7 = 1 (right panel). 
The dashed curves corresponds to the power law t~ 2//3 for 7 = N and to t^ 1 for 
7 = 1. 



positive and negative velocity seems to move ballistically. On the other hand, 
when 7 = N, the group velocity of the outgoing peaks depends weakly on time, 
probably reaching a final constant velocity at much longer times. 

The scalings (6.2) and (6.3) are valid for the decaying of the initial perturba- 
tion, namely they are valid if measured from the decay of the central peak. Nev- 
ertheless, we find that in our case, similar scalings are possible for the outgoing 
peaks. Assuming that the excess of energy is transported across the system as a 
density packet, whose area is preserved on average, we show in Fig. 10 the damp- 
ing of the amplitude of the moving peak as a function of time. For 7 = N (left 
panel), the amplitude of the peak decays as t -2 / 3 , corresponding to a heat con- 
ductivity that scales as k ~ iV 1 / 2 . Note, however, that within numerical accuracy, 
the decay could also be consistent with k ~ iV 1 / 3 as it is found at the beginning 
of this section. In fact, from a fit to a power law of the amplitude decay we have 
obtained scalings for the amplitude decay between 0.62 to 0.66. In any case, the 
spreading of the outgoing peak reflects, locally, the anomalous character of the 
heat transport. On the other hand, the amplitude decay for 7 = 1 (right panel) is 
consistent with t~ l , corresponding to a heat conductivity that grows linearly with 
N, in full agreement with the observations at the beginning of this section. 
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